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Abstract. We present quasi-stationary simulations of three-dimensional models with 
a single absorbing configuration, viz. the contact process (CP), the susceptible- 
infected-susceptible (SIS) and the contact replication process (CRP). The moment 
ratios of the order parameters for DP class in three dimensions were set up using the 
well established SIS and CP models. We also show that the mean-field exponents in 
d = 3 reported previously for CRP [Ferreira SC 2005 Phys. Rev. E 71 017104] is a 
transient observed in the spreading analysis. 
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1. Introduction 

Phase transitions to a single absorbing configuration, a state in which the system can not 
scape from, are nowadays a topic in the frontier of Nonequihbrium Statistical Physics 
[HE]. Concomitantly with the increasing interest on absorbing/active phase transitions 
in complex topologies [31 lU El El [7] , there are still a lot of open problems being intensively 
investigated on regular lattices such as the effects of quenched disorder [H El [10] , diffusion 

as well the modeling of predator-prey systems [I2], and clonal replication [T31 ITi] . 

Under the renormalization group point of view, it is expected [H [151 US] that the 
absorbing phase transitions in models with a positive one-component order parameter, 
short-range interactions and without additional symmetries or quenched disorder belong 
generally to the universality class of directed percolation (DP). This conjecture is known 
as Janssen-Grassberger criterion [1]. It is worthwhile to mention, the interest on this 
kind of phase transitions was raised by the recent experimental observation of the DP 
class in absorbing-state phase transitions [HI |19]. On the other hand, while DP is 
considered the most robust universality class of the absorbing-state phase transitions, 
the precise numerical determination of the critical exponents of a specific model can be 
masked by factors like diffusion [TT] and weak quenched disorder p!0]- 

The contact process (CP), the standard example of the DP universality class, is a 
toy model of epidemics [20] More recently, a novel variation of the CP was introduced 
for the modeling of clonal (copies of themselves) replication, the contact replication 
process (CRP) [131 [Tl]. Since neither additional symmetries nor long-range interactions 
were included, the CRP fulfils the requirements of the Janssen-Gassberger criterion. 
However, the first dynamic spreading analysis of CRP in ci = 1 — 3 dimensions, reported 
in [131 E] intriguingly classified the model in the DP universality class in one and two, 
but not in three dimensions. Surprisingly, in d = 3 the reported spreading exponents 
were those predicted by the mean-field approach [H]. 

In the present work we applied spreading analysis and the method of quasi- 
stationary simulations [22l [23] in three-dimensional models that fulfill the Janssen- 
Grassberger criterion. Particularly, we turned back to the CRP model and showed 
that the mean-field behaviour observed previously in d = 3 [H] is a transient associated 
to the closeness between critical creation and annihilation events. Additionally, we have 
analysed several moment ratios in d = 3 and have determined the universal values for 
DP class based on the results obtained from CP and SIS models. The paper is outlined 
as follows. In section [21 the models and simulation procedures are described. Simulation 
results are presented and discussed in section [3l Conclusions are drawn in section [H 

X A CP variation called susceptible-infected-susceptible (SIS) model is more widely applied in 
epidemiological studies |3]. 
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2. Models and methods 

2.1. Models 

The contact process (CP) is defined in a hypercubic lattice Z'^ in which the sites represent 
individuals in two states: healthy {Zi = 0) or infected {Zi = 1). Infected sites become 
healthy (empty) at unitary rate 1 while healthy sites are infected at a rate UiX/q, where 
q is the lattice coordination number and rii the number of infected sites surrounding 
(first-neighbours in the distance) the healthy site i. These rules implies, for values of 
infection rate above a certain Ac, an infection fiowing equally distributed among all 
nearest neighbours (NN) of infected sites. 

The susceptible-infected- susceptible model is a variation of the CP dynamics in 
which any empty site with one or more infected nearest-neighbours becomes infected at 
rate A. Again, a unitary cure rate is assumed. In the literature, the SIS model is also 
known as the A model |17j . 

Contact replication process rules are are very similar to those of CP. Instead of 
individuals, the sites represent places where cells lie. Analogously to the spontaneous 
cure in CP, cells die at unitary rate. However, a cell replicates at a rate A and the 
offspring occupies one of its empty NN chosen at random. So, an empty site i is occupied 
at rate A Zj/uj. The sum is done over all neighbours of the site i and Uj is defined 
as before. Notice that the creation process is facilitated in the CRP in relation to 
CP since the occupation fiows uniformly among only empty neighbours implying lower 
critical rates. Indeed, estimates reported for CRP were Ac = 2.02634(4), 1.08320(7), 
and 1.0000(1) for d = 1,2, and 3 fi3\ HI], in comparison with those for the CP 
Ac = 3.29785(2), 1.64877(3), and 1.31686(1) pi], respectively. Notice that the three 
models share the same symmetries and, consequently, they are expected to belong to 
the same universality class. 

For all these models, Monte Carlo simulations were performed using the usual 
procedure [1]: An event, creation or annihilation, is selected with probabilities p = 
A/(l + A) and I — p, respectively, and an occupied site i is chosen at random. In the 
annihilation process, the occupied site become empty in all models while the creation 
depends on model. In CP, one nearest-neighbour (NN) of i is chosen at random and 
infected if empty, otherwise nothing occurs. In SIS, all empty sites neighbouring the site 
i are occupied. Finally, in CRP one of empty neighbours, if there are anyone, is chosen 
at random and occupied. In all cases, the time is incremented by At = 1/A^, where N 
is total number of infected sites. 

2.2. Quasi-stationary simulations 

Stationary analysis of systems with transitions to absorbing configurations in the 
proximity of the critical point are ruled by strong finite size effects. Indeed, the unique 
actual stationary state of finite systems is the absorbing one. A common alternative to 
avoid this difficulty is to restrict the averages to the survival samples and apply a finite 
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size analysis. However, such procedure is not free of ambiguities or misinterpretations 
|23j . An alternative approach is the quasi-stationary QS simulation method [22] 123] . 
This method consists of storing a list with M configurations visited in the history of 
the system and periodically replacing one of them by the current state. Whenever 
the system try to visit the absorbing state, the configuration is replaced by an active 
one selected at random from the list containing the sample of configurations and the 
simulation continues as usually. 

The QS simulations were performed as follows. Firstly, the list of configurations is 
incremented whenever the time increases by a unity up to a list with M configurations 
is achieved. Secondly, a configuration of the list randomly chosen is replaced by the 
current one with a given probability prep- We used a large value of Prep = 0.05 for a 
initial relaxation period, more precisely for t < 5 x 10^ aiming to speed up the erasing 
of the memory of the initial conditions. In turn, prep = 2 x 10~^ was adopted for 
the remaining of the simulation. Runs with = 2 x 10*^ steps and averages after a 
relaxation time tr = 1 x 10^ were used. Finally, the averages and uncertainties were 
obtained with at least 15 (to the largest system) independent runs with a full lattice as 
initial condition. Periodic boundary condition were always used. 

2.3. Spreading analysis 

The critical point Ac can be efficiently determined by the spreading analysis, which 
consists of evolving the system from a perturbation to the absorbing state (a single 
occupied site at the origin of the lattice) and computing the survival probability P and 
mean number of occupied sites as time functions. In this analysis, the averages are 
done over all samples, surviving or not. Asymptotic power law dependencies, 

P{t) ~ t-^ and N{t) ~ f (1) 

are expected at criticality. Since N at upper-critical and exponentially decays in 
the sub-critical regimes, deviations from power laws are expected around the critical 
point and a null curvature criterion of the double-logarithm plots of A^ against t can 
be used to determine the critical point [13]. Analogous analysis can be done with the 
survival probability. 

The spreading is defined by 




in which rj is the distance from the original seed where the perturbation was introduced. 
At criticality, it is expected an asymptotic power law 

R\t) ~ t'. (3) 

The spreading exponents obey a hyperscahng relation 46 + 2ri = dz [H [2] . 

In order to illustrate the method, we present results from quasi-stationary 
simulations of the two-dimensional CRP model§|. The QS density vanishes as usual 

§ in rcf. [14^ only results from spreading simulations were reported. 
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Figure 1. Critical QS density ps and lifetime t versus system size L for the two- 
dimensional CRP model. Inset: critical moment ratio m = (p^)/(p)^ for varying 
system sizes. Dashed line is the moment ratio for the contact process taken from [17) . 



at A = Ac = 1.08322(2) when L — t- cxo (the critical rate were determined using the 
spreading analysis). Figure [1] shows the critical densities and hfetimes of the QS state 
versus system size. Lifetime was calculated as [22j 

r = IM, (4) 

where pi is the probability of attempting the absorbing configuration in QS simulations. 
The power laws 

Ps ~ L-^/-"^ (5) 

and 

r ~ L'^li/'^^ (6) 

were verified and the slopes for L > 20 provide = 0.800(7) and i^\\/i^± = 1.764(14). 
As expected, these exponents are in very good accordance with those reported for the 
DP universality class [21] (tabled]). Also, the critical variance of the order parameter, 
defined by x = L^iip"^) — (p)^) [Ij? scales as 

X-^L-r/^^, (7) 

with an exponent 7/z/_l = 0.412(8), again agreing with the DP value. 

The ratio between moments of the order parameter is widespread as an useful tool 
for the determination of the critical point of equilibrium systems [21]. This concept 
has been extended to non-equilibrium systems, particularly to models belonging to the 
DP universality class for which the moment and cumulant ratios have proven to be 
universal quantities [IT]. The inset of figure IH^b) shows the ratio m = (p^)/(p)^ for 
the two-dimensional CRP model at criticality for varying system sizes. It was shown 
that this ratio assumes the critical value nic = 1.3257(5) for CP in two dimensions [IT] . 
Using data for L > 80 we found nic = 1.3274(9) for two-dimensional CRP, again in 
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accordance with DP universality class. Additionally, we also verified the agreement of 
higher moment ratios with the predictions for DP classes for d = 1 and 2. 

Table 1. Critical exponents for DP class calculated from spreading exponents given in 
[2T . Scaling relations /? — Sv^, z — ^v^jv^^ and 7 = dv\_ — 2^ [11 were used whenever 
necessary. 



Quasistationary 


Spreading 




Exponent d = 2 d = 3 


Exponent d = 2 


d = 3 


u\\/vi_ 1.765(3) 1.919(4) 


S 0.452(1) 


0.756(1) 


fi/vi_ 0.799(2) 1.394(5) 


77 0.229(3) 


0.110(1) 


-l/vx_ 0.401(4) 0.212(96) 


z 1.133(2) 


1.042(2) 




Figure 2. Left: critical QS quantities against system size for three-dimensional CP. 
Lines are least square fits by ([5]). Right: critical QS quantities for SIS rescaled by pure 
(filled symbols) and corrected (open symbols) power laws. Curves were shifted for sake 
of visibility. 

Quasistationary simulations of the CP and SIS in d = 3 dimensions are shown 
in figure [2l The CP critical rate Ac = 1.31686(1) was taken from [21] whereas 
Ac = 0.24805(2) for the SIS model was estimated in the present work using spreading 
analysis. Least square fits for L > 20 provide the exponents (3/h'± = 1.400(8), 7/i^± = 
0.233(5), and z/y/i^x = 1-919(9) for CP while = 1.395(5), = 0.252(14), and 

'^\\/'^± = 1.944(6) were obtained for SIS model. Notice that all of them are consistent 
with DP class, although data are deviated from pure power laws for small sizes. This 
effect can be taken in to account with a suitable correction to the amplitude of the scale 
law for lifetime given by [23j 

Ml A 

In r = — In L + — ^ + const., (8) 
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and so on for the others QS quantities. Actually, our results are not significantly affected 
by the particular choice of the correction. The value = 0.75 was established for CP 
in d = 1 [23] and also adopted here for = 3. Using this correction, the exponents 
are /3/i'± = 1.395(4), = 0.216(3), and i^\\/i^± = 1.916(5) for contact process 

representing a significant reduction of the error estimates and an improvement of the 
closeness with DP class. For the SIS model, the correction to the scaling provides 
l3/ui_ = 1.403(4), 7/z/_L = 0.209(2), and z/||/z/^ = 1.922(2), an even better improvement 
of the estimates as well as their proximities with DP class. In left panel of figure El the 
rescalings by pure power laws, t/L'^^^^'^^ and so on, are compared with those given by 
dH]). It is neat the overmatching fit between data and the ansatz (jH]). 




10' 10^ 10' 10* 

t 



Figure 3. Spreading analysis for the tlrree-dimensional CRP around tire critical point. 
Bottom and top groups of curves correspond to survival probability P and mean 
number of occupied sites N. In each group, are shown curves for A — 1.00362, 1.00363, 
and 1.00364 from bottom to top. Dashed lines represent slopes and -1 corresponding 
to the DP mean field exponents 77 and S, respectively. Inset: mean critical density 
(not the QS one) rescaled by the power law for a initial condition with all sites 
occupied. 

A central point of the present work is the QS analysis of CRP in three-dimensions 
for which the universality class was left as an open question. Firstly, we recalculated the 
critical point determination using spreading analysis and found a value approximately 
0.3% larger than the originally reported pH]. Using the criterion of downward and 
upward curvatures of the plots N vs. t, our best estimate is Ac = 1.00363(1). This small 
discrepancy is central since this certainly excludes the equality between creation and 
annihilation critical rates. The curves N{t) and P{t) have transients consistent with 
the DP mean-field exponents (figure [3]), in agreement with the previously reported CRP 
simulations [H]. Even though a large power interval was not obtained for the largest 
time simulated, the exponents fitted in the interval t = 10^ — 10^, namely, rj = 0.09(2), 
S = 0.78(3) and z = 1.04(1), do not exclude the DP class (table [1]). Additionally, 
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the critical density for a fully occupied initial configuration decays as p ~ t o.76(i) 
t > 10^, in excellent agreement with DP class (inset of figure [3]). 



10 
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Figure 4. QS simulations of CRP model at criticality. Solid lines are least square fits 
using ([5|) and dashed lines are power laws with the DP exponents taken from table [TJ 



Results of QS simulations for the CRP model at the critical point are shown in 
Fig. HJ The QS exponents obtained from fits in the range L > 20 are = 1.325(6), 
7/z/_L = 0.29(5), and i'\\/i'± = 1.79(1) which are not concluding about the agreement 
with the DP universality class. However, with exception of the lifetime, the exponents 
seem converging to the DP value. Applying a correction to the scaling (|8]), the corrected 
exponents for ps and X : f^/^± = 1-374(3) and 7/z/_l = 0.19(4), are closer to DP class. 
The lifetime analysis deserves some comments. The QS method used in present work, 
involves the transition to the absorbing configuration passing by the state with a single 
occupied site. In the critical CRP, the creation event occurs with a frequency slightly 
larger than the annihilation, which can be easily verified from the model critical rate 
together with its rules. Consequently, large clusters of occupied sites are much more 
seldom in CRP than in CP and, mainly, in SIS models implying in a too small frequency 
of visiting the pre-absorbing state (figure |5]). Moreover, the average time that a sample 
is kept in the list containing the system history, given by M/prep, must be much larger 
than the sample lifetime, r = 1/pi, to avoid a same run anomalously contributing many 
times to the list of configurations [23]. Obviously, this effect is enhanced as the system 
size increases. Thus, a computationally prohibitive small Prep and/or large M together 
with too large relaxation and averaging times are demanded. Instead of lifetime, we can 
determine the characteristic time Tp to the density reaches the QS state in conventional 
QS simulations [H [13]. The scale relations p ~ L^-^'^C^) and ~ were found for 

CRP, in agreement with the DP class as predicted by Janssen-Grassberger conjecture. 

Moments and cumulants for DP class in d = 1 — 5 with a homogeneous constant 
external source were recently reported by Janssen et al. [25j. So, let us introduce the 



Quasi- stationary simulations of DP class in d = 3 



9 




3.0 



K 2.0 



1.0 



O-O SIS 
CB-El CP 
CRP 




G — o ooooooooooooooooooooooo 



10 



10 



10 

t 



10 



10 



Figure 5. Critical quantities normalized by the SIS model. Top: lifetimes as functions 
of the system sizes. Bottom: squared spreading radius versus time. As one can see, 
diffusion and lifetimes are enhanced in CRP model. 



notation yU„ = (p") and /€„ to the nth moment and cumulant of the order parameter, 
respectively. In particular, we have 

K2 = fl2- (A (9) 

and 

K4 = /i4 — 4yU3yUi — 3/^2 + 12/U2Aii — 6p^. (10) 

Dickman and da Silva [17J showed that the ratios ^2/ /^s//^?, /^3//^i/^25 and 
fi4/fj^2 assume universal values for DP class in d = 1 and d = 2 dimensions. It is 
important to notice that both even and odd moments can be studied, since the order 
parameter is non- negative. Several moment and/or cumulant ratios as functions of the 
system size for 3d models are shown in figure El The moment ratios for SIS and CP 
models exhibit finite size effects, but seem to monotonically converge to a constant 
value when L — )■ 00. Assuming a correction given by m(L) = m(oo) + aL~^, where ip 
and a are fit parameters, the asymptotic moment ratios can be estimated. Since some 
data for fourth moment were ruled by very large error bars, these data were excluded 
from the fits. Extrapolated moment ratios for CP and SIS are listed in table [2] and are 
consistent with the hypothesis of universality. It is worth to stress that ratios shown in 
table [2] differ from those reported by Janssen et al. |25] since, in our studies, there is no 
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Figure 6. Moment rations for SIS, CP and CRP models in d = 3 dimensions. Solid 
lines are non-linear fits by function m{L) = 171(00) + aL~'^ with a correlation coefficient 
r > 0.999. 



external source. The same occurs in d = 1 and 2 when the results of Janssen et al. 
are compared with those from reference |17j . 

CRP seems to be different from CP and SIS. Indeed, instead of the monotonic 
convergence to a constant value, the moment ratios decreases after a maximum. This 
difference can be again associated to the proximity between critical creation and 
annihilation events and the resultant high diffusivity. Actually, the ratio (^2/^^! = 
1^2 1 + 1 first grows towards the value 1.550. This value is nearby the established 
value of 1.660 for CP on a complete graph [22]. Again, we have a mean- field behaviour 
for small system that must converge to the usual DP for asymptotic large systems. 



Table 2. Asymptotic moment ratios for SIS and CP in c? = 3. 



Model 


K2/m\ 










SIS 
CP 


0.469(3) 
0.470(2) 


0.490(6) 
0.454(10) 


2.678(12) 
2.697(6) 


1.822(3) 
1.833(3) 


2.629(12) 
2.649(18) 



4. Conclusions 

We performed large-scale simulations of the contact process, the susceptible-infected- 
susceptible and the contact replication process in three dimensions. Applying the quasi- 
stationary simulation method, we were able to determine the moment ratios of the order 



Quasi- stationary simulations of DP class in d = 3 



11 



parameters for DP class in three dimensions in the absence of an external field. We also 
show that the mean-field exponents in c? = 3 for the CRP, reported in p!^ are a transient 
observed in the spreading analysis. 

Quasi-stationary simulations and suitable corrections to the scaling revealed that 
the CRP model belongs to the directed percolation universality class, as expected by 
the Janssen-Grassberger criterion. However, the moment ratios for CRP agree with the 
universal DP values in ci = 1 and 2 dimensions but do not in d = 3. The discrepancy 
lies on the closeness between critical creation and annihilation events, which gives rise 
to a diffusive transient behaviour in the CRP model, also responsible by the transient 
mean field exponents obtained in ci = 3. 
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